#####################################################################################
# Models for Main Results

m.comp.1 <- lm(liberal ~ rel_banperm + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.1)
y.comp.1 <- lm(rel_feel ~ rel_banperm*liberal + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.1)

m.comp.2 <- lm(liberal ~ rel_banperm + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.2)
y.comp.2 <- lm(rel_feel ~ rel_banperm*liberal + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.2)

m.comp.3 <- lm(liberal ~ rel_banperm + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.3)
y.comp.3 <- lm(rel_feel ~ rel_banperm*liberal + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.3)

m.comp.4 <- lm(liberal ~ rel_banperm + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.4)
y.comp.4 <- lm(rel_feel ~ rel_banperm*liberal + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.4)

m.comp.5 <- lm(liberal ~ rel_banperm + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.5)
y.comp.5 <- lm(rel_feel ~ rel_banperm*liberal + fem + I(age/10) + edu_high + religiosity + leftright + rel_muslim + rel_dev + rel_rad + rel_immigrant, data=imp.data.5)

c.1 <- coef(m.comp.1)
c.2 <- coef(m.comp.2)
c.3 <- coef(m.comp.3)
c.4 <- coef(m.comp.4)
c.5 <- coef(m.comp.5)

s.1 <- se.coef(m.comp.1)
s.2 <- se.coef(m.comp.2)
s.3 <- se.coef(m.comp.3)
s.4 <- se.coef(m.comp.4)
s.5 <- se.coef(m.comp.5)

coefs <- cbind(c.1, c.2, c.3, c.4, c.5)
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(s.1^2, s.2^2, s.3^2, s.4^2, s.5^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

s.combined <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, s.combined)
m.combined.res <- round(combined.res, 2)

m.combined.res <- rbind(m.combined.res, display(m.comp.1)$n, display(m.comp.1)$r.squared)
m.combined.res <- round(m.combined.res, 2)

colnames(m.combined.res) <- c("est", "se")
rownames(m.combined.res)[2] <- c("liberal_policy")
rownames(m.combined.res)[12] <- c("N")
rownames(m.combined.res)[13] <- c("R2")

c.1 <- coef(y.comp.1)
c.2 <- coef(y.comp.2)
c.3 <- coef(y.comp.3)
c.4 <- coef(y.comp.4)
c.5 <- coef(y.comp.5)

s.1 <- se.coef(y.comp.1)
s.2 <- se.coef(y.comp.2)
s.3 <- se.coef(y.comp.3)
s.4 <- se.coef(y.comp.4)
s.5 <- se.coef(y.comp.5)

coefs <- cbind(c.1, c.2, c.3, c.4, c.5)
c.combined <- apply(coefs, 1, mean)

var.within <- cbind(s.1^2, s.2^2, s.3^2, s.4^2, s.5^2)
var.within <- apply(var.within, 1, mean)

m <- 5
var.between <- (coefs - c.combined)^2
var.between <-  apply(var.between, 1, sum)
var.between <-  (m-1)^-1 * var.between

s.combined <- sqrt(var.within + (1 + m^-1) * var.between )

combined.res <- cbind(c.combined, s.combined)
y.combined.res <- round(combined.res, 2)


y.combined.res <- rbind(y.combined.res, display(y.comp.1)$n, display(y.comp.1)$r.squared)
y.combined.res <- round(y.combined.res, 2)

colnames(y.combined.res) <- c("est", "se")
rownames(y.combined.res)[2] <- c("liberal_policy")
rownames(y.combined.res)[3] <- c("citizen_reaction")
rownames(y.combined.res)[13] <- c("liberal_policy*citizen_reaction")
rownames(y.combined.res)[14] <- c("N")
rownames(y.combined.res)[15] <- c("R2")
